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Abstract:  This  paper  deals  with  the  formulation,  calibration,  and  validation  of  the  Lattice  Discrete  Particle 
Model  (LDPM)  suitable  for  the  simulation  of  the  failure  behavior  of  concrete.  LDPM  simulates  concrete  at  the 
meso-scale  considered  to  be  the  length  scale  of  coarse  aggregate  pieces.  LDPM  is  formulated  in  the  framework 
of  discrete  models  for  which  the  unknown  displacement  field  is  not  continuous  but  only  defined  at  discrete  points 
representing  the  center  of  discrete  particles.  The  size  and  distribution  of  the  particles  is  obtained  by  idealizing 
the  geometry  of  the  concrete’s  internal  structure.  Discrete  compatibility  and  equilibrium  equations  are  used  to 
formulate  the  governing  equations  of  the  LDPM  computational  framework.  Particle  contact  behavior  represents 
the  mechanical  interaction  among  adjacent  aggregate  particles  through  the  embedding  mortar.  Such  interaction 
is  governed  by  meso-scale  constitutive  equations  simulating  meso-scale  tensile  fracturing  with  strain-softening, 
cohesive  and  frictional  shearing,  and  nonlinear  compressive  behavior  with  strain-hardening.  The  present.  Part  I, 
of  this  two-part  study  deals  with  model  formulation  leaving  model  calibration  and  validation  to  the  subsequent 
Part  11. 

1  Introduction 

Concrete  is  an  heterogeneous  material  characterized  by  several  length  scales  of  observation  ranging  from 
the  atomistic  scale  m),  characterized  by  the  behavior  of  crystalline  particles  of  hydrated  Portland 

cement,  to  the  macroscopic  scale  (10^  m),  at  which  concrete  has  been  traditionally  considered  homoge¬ 
neous.  It  is  now  widely  recognized  that  accurate  modeling  of  multiscale  materials  calls  for  the  adoption 
of  multiscale  techniques  able  to  bridge  the  various  scales  and  to  bring  to  the  macroscopic  scale  the  most 
important  effects  of  lower  scale  phenomena.  In  the  recent  past,  publications  proposing  new  multiscale 
theories  have  flourished,  especially  for  modeling  nano-composite  materials  and  atomistic  and  molecular 
systems  [23].  The  same  kind  of  development  has  not  appeared  yet  in  concrete  mechanics  literature  and 
in  civil  engineering  in  general.  The  main  reason  for  this  can  be  traced  back  to  the  extreme  complexity  of 
concrete  internal  structure  and  to  the  unavailability  of  accurate  hne-scale  models  for  concrete. 

In  the  last  twenty  years,  various  authors  attempted  the  development  of  concrete  models  targeting 
concrete  mini-scale  (length  scale  of  10“‘^m  or  less)  and  meso-scale  (length  scale  10“^m).  The  mini-scale 
nomenclature  was  first  introduced  by  Cusatis  et  ah  [19]  and  is  relevant  to  the  description  of  concrete  as  a 
three-phase  material:  cement  paste,  aggregate,  and  interfacial  transitional  zone,  whereas  the  meso-scale  is 


2 


relevant  to  the  characterization  of  concrete  as  two-phase  material:  mortar  and  coarse  aggregate.  It  must 
be  noted  that  some  authors  use  the  term  “meso-scale”  in  a  wider  sense  to  include  the  “mini-scale” . 

Mini-scale  models  were  proposed  by  several  authors  [31,  33,  32,  10,  1,  9,  37].  Remarkable  are  the 
contributions  due  to  Wittmann  and  coworkers  [31]  for  2D  models,  and  to  Carol  and  coworkers  [12,  11,  13] 
for  3D  models.  They  used  hnite  element  techniques  to  model,  with  different  constitutive  laws,  coarse 
aggregate  pieces,  mortar  matrix,  and  an  inclusion-matrix  interface.  A  more  effective  alternative  to  the  use 
of  hnite  elements  was  proposed  by  Van  Mier  and  coworkers  [33,  32]  who  removed  the  continuum  hypothesis 
and  modeled  concrete  through  a  discrete  system  of  beams  (lattice).  In  their  approach,  lattice  meshes  were 
superimposed  to  digitalized  images  of  the  concrete  internal  structure  to  assign  different  material  properties 
to  the  lattice  elements  corresponding  to  the  various  components  (matrix,  aggregate,  and  interface).  Along 
this  line,  Bolander  and  coworkers  [9,  37]  formulated  a  discrete  mini-scale  model  based  on  the  interaction 
between  rigid  polyhedral  particles  obtained  though  the  Voronoi  tessellation  of  the  domain.  Mini-scale 
models  provide  realistic  simulations  of  concrete  cracking,  coalescence  of  multiple  distribute  cracks  into 
localized  cracks,  and  fracture  propagation.  However,  they  tend  to  be  computationally  intensive  especially 
for  3D  modeling  that  is  required  to  correctly  capture  compressive  failure  and  conhnement  effects. 

Computationally  less  demanding  are  the  meso-scale  models  [5,  17,  19]  in  which  the  basic  material 
components,  whole  aggregate  pieces  and  the  layer  of  mortar  matrix  between  them,  are  modeled  through 
discrete  elements  (either  lattice  elements  or  discrete  particles)  but  are  themselves  not  discretized  on  a 
hner  scale.  Meso-scale  models  greatly  reduce  the  size  of  the  numerical  problems  but  at  the  same  time 
can  capture  the  fundamental  aspects  of  material  heterogeneity.  Meso-scale  models  have  made  possible 
the  realistic  simulation  of  both  tensile  and  compressive  softening.  Preliminary  results  on  the  modeling  of 
multiaxial  behavior  and  conhnement  effects  were  also  achieved  by  Cusatis  et  ah  [18]  and  Belheine  et  ah  [7]. 
The  main  objective  of  this  article  is  to  discuss  a  recently  developed  meso-scale  model  for  concrete,  called 
the  Lattice  Discrete  Particle  Model  (LDPM).  The  development  of  LDPM  is  a  synthesis  of  two  independent 
research  efforts  that  led  to  the  formulation  of  the  Conhnement  Shear  Lattice  (CSL)  Model  [19,  17,  18]  and 
the  Discrete  Particle  Model  (DPM)  [26]. 

LDPM  shares  the  following  features  with  CSL;  (a)  It  simulates  concrete  mesostructure  by  a  system 
of  interacting  aggregate  particles  connected  by  a  lattice  system  that  is  obtained  through  a  Delaunay 
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tetrahedralization  of  the  aggregate  centers;  (b)  The  position  of  each  aggregate  piece  throughout  a  given 
concrete  specimen  is  dehned  by  means  of  the  basic  concrete  properties  and  the  size  distribution  of  the 
aggregates;  (c)  The  geometrical  interaction  between  the  particles  is  obtained  by  a  three-dimensional  domain 
tessellation  defining  a  set  of  polyhedral  cells  each  including  one  aggregate  piece;  (d)  The  mechanical 
interaction  between  the  particles  is  characterized  by  both  normal  and  shear  stresses;  and  (e)  The  meso- 
scale  constitutive  behavior  is  softening  for  pure  tension  and  shear-tension  while  it  is  plastic  hardening  for 
pure  compression  and  shear-compression. 

LDPM  inherited  from  DPM  the  MARS  (Modeling  and  Analysis  of  the  Response  of  Structures)  com¬ 
putational  environment  [27,  28,  29]  that  includes  long  range  contact  capabilities  typical  of  the  classical 
formulation  of  Discrete  Element  Methods  (DEM)  [15].  This  feature  is  particularly  important  for  simulating 
pervasive  failure  and  fragmentation. 

While  building  on  the  successful  developments  of  CSL  and  DPM,  LDPM  formulation  is  characterized 
by  a  number  of  new  features  that  greatly  enhance  its  modeling  and  predictive  capabilities.  These  new 
features  can  be  summarized  as  follows: 

1.  Interaction  among  the  particles  is  formulated  through  the  analysis  of  an  assemblage  of  four  aggregate 
pieces  whose  centers  are  the  vertexes  of  the  Delaunay  tetrahedralization.  This  makes  possible  the 
inclusion  of  volumetric  effects  in  the  constitutive  law  that  cannot  be  taken  into  account  by  the 
two-particle  interaction  used  in  CSL  and  DPM. 

2.  Stresses  and  strains  are  defined  at  each  single  facet  of  the  polyhedral  cells  containing  the  aggregate 
pieces.  Compared  to  previous  formulations,  this  allows  a  better  stress  resolution  in  the  mesostructure, 
which,  in  turn,  leads  to  a  better  representation  of  meso-scale  crack  and  damage  distribution. 

3.  The  constitutive  law  simulates  the  most  relevant  physical  phenomena  governing  concrete  damage 
and  failure  under  tension  as  well  as  compression.  Compared  to  the  constitutive  law  used  in  the 
previous  work  [17],  the  present  law  provides  better  modeling  and  predictive  capabilities  especially 
for  the  macroscopic  behavior  in  compression  with  conhnement  effects. 

4.  The  constitutive  equations  include  simple  but  effective  unloading-reloading  rules  that  permit  an 
accurate  simulation  of  concrete  response  under  cyclic  loadings. 
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5.  The  constitutive  equations  also  include  the  effect  of  material  compaction  and  densification  due  to 
the  effect  of  high  confining  pressures. 

The  present  formulation  can  realistically  simulate  all  aspects  of  concrete  response  under  quasi-static 
loading,  including  tensile  and  compressive  strength,  cohesive  fracture  and  size  effect,  damage  in  compres¬ 
sion,  compression-shear  behavior  with  softening  at  zero  or  low-confinement  and  hardening  at  high  confined 
compression,  and  strength  increase  under  biaxial  loading.  This  paper  (Part  I  of  a  two-part  study)  discusses 
the  details  of  LDPM  formulation  and  its  numerical  implementation.  Part  II  will  focus  on  its  extensive 
calibration  and  validation. 

2  Geometrical  Characterization  of  Concrete  Mesostructnre 

The  geometrical  characterization  of  concrete  mesostructnre  is  based  on  a  four-step  procedure  that  aims  at 
defining  (1)  the  number  and  size  of  coarse  aggregate  pieces  (particles);  (2)  particle  position;  (3)  interparticle 
connections;  and  (4)  surfaces  through  which  forces  are  transmitted  between  adjacent  particles.  These 
surfaces  will  also  represent  weak  locations  in  the  concrete  meso-structure  where  damage  is  likely  to  localize. 

2.1  Particle  Generation 

In  the  first  step,  particle  generation  is  carried  out  by  assuming  that  each  aggregate  piece  can  be  approxi¬ 
mated  as  a  sphere.  Under  this  assumption,  typical  concrete  granulometric  distributions  can  be  represented 
by  the  particle-size  distribution  function  (psd)  proposed  by  Stroeven  [35]: 

=  [1  -  {do/l)^]d^+^ 

where  da  is  the  maximum  aggregate  size,  and  do  is  the  minimum  particle  size  used  in  the  simulations, 
do  7^  0  to  limit  the  number  of  degrees  of  freedom.  The  above  psd  can  be  interpreted  as  the  probability 
density  function  (pdf)  for  the  occurrence  of  a  certain  diameter  d.  The  cumulative  distribution  function 
(cdf)  can  be  then  computed  as 


P(d) 


/(d)dd 


'do 


1  ~  {do/ dy 
1  -  {do/daY 


(2) 
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It  can  be  shown  [35]  that  the  psd  in  Equation  1  is  associated  with  a  sieve  curve  (percentage  of  aggregate 
by  weight  retained  by  a  sieve  of  characteristic  size  d)  in  the  form 

where  np  =  3  —  q-  For  q  =  2.5  {np  =  0.5),  Equation  3  corresponds  to  the  the  classical  Fuller  curve  which 
for  its  optimal  packing  properties,  was  extensively  used  in  concrete  technology  [24]. 

For  a  given  cement  content  c,  water-to-cement  ratio  w/c,  specimen  volume  V,  maximum  aggregate  size 
da,  and  minimum  particle  size  do  (which  governs  the  resolution  of  the  model),  particles  to  be  placed  inside 
the  volume  can  be  obtained  as  follows: 

1.  Compute  aggregate  volume  fraction  as  =  1  —  cj pc  —  w/ —  Vair,  where  w  =  {w/c)c  is  the  water 
mass  content  per  unit  volume  of  concrete,  w/c  is  the  water-to-cement  ratio,  pc  =  3150  kg/m^  is 
the  mass  density  of  cement,  p^  =  1000  kg/m^  is  the  mass  density  of  water,  and  Vair  is  the  volume 
fraction  of  entrapped  or  entrained  air  (typically  3-4%); 

2.  Compute  the  volume  fraction  of  simulated  aggregate  as  Vao  =  [1  —  F{dQ)\va  =  [1  —  {d^/ daY^]va', 

3.  Compute  the  total  volume  of  simulated  aggregate  as  Vao  =  VaoV; 

4.  Compute  particle  diameters  by  sampling  the  cdf  in  Equation  2  by  a  random  number  generator: 
di  =  do[l  —  Pj(l  —  dg/d^)]”^/'^,  where  Pi  is  a  sequence  of  random  numbers  between  0  and  1.  Figure 
la  shows  a  graphical  representation  of  the  particle  diameter  selection  procedure. 

5.  For  each  newly  generated  particle  in  the  sequence,  check  that  the  total  volume  of  generated  particles 
%ao  =  X^i('^'^?/6)  does  not  exceed  I4o-  When,  for  the  first  time,  Wo  >  Wo  occurs,  the  current 
generated  particle  is  discarded,  and  the  particle  generation  is  stopped. 

Figure  lb  shows  the  comparison  between  the  theoretical  sieve  curve  (solid  line)  and  the  computational 
sieve  curve  (circles),  obtained  through  the  procedure  highlighted  above  for  the  generation  of  a  100-mm-side 
cube  of  concrete  characterized  by  c  =  300  kg/m^,  w/c  =  0.5,  np  =  0.5,  do  =  4  mm,  and  da  =  S  mm. 

In  order  to  simulate  the  external  surfaces  of  the  specimen  volume,  the  generated  particles  are  augmented 
with  zero-diameter  particles  (nodes).  Assuming  that  the  external  surfaces  of  the  specimen  volume  can  be 
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described  through  a  set  of  vertexes,  edges,  and  polyhedral  faces,  one  node  for  each  vertex  is  hrst  added 
to  the  particle  list.  Then,  iVe  =  INT(Le/h5)  and  Np  =  INT(Ap/hf)  (where  INT  is  the  integer  part  of  the 
number)  nodes  are  associated  with  each  edge  e  and  polyhedral  face  p,  respectively.  Le  is  the  length  of  a 
generic  surface  edge,  Ap  is  the  area  of  a  generic  surface  polyhedron,  and  the  average  surface  mesh  size  hg 
is  chosen  such  that  the  resolution  of  the  discretization  on  the  surface  is  comparable  to  the  one  inside  the 
specimen.  This  is  achieved  by  setting  hg  =  ^gdo.  Numerical  experiments  performed  in  this  study  show 
that  ^g  =  1.5  provides  optimal  results. 

2.2  Definition  of  Particle  Position 

The  second  step  of  the  geometrical  characterization  of  concrete  mesostructure  is  relevant  to  the  random 
distribution  of  the  generated  particles  on  vertexes,  edges,  surface  faces,  and  interior  volume.  First,  the 
vertex  nodes  are  placed.  Secondly,  nodes  are  distributed  over  the  edges  and  surfaces  by  allowing  a  minimum 
distance  of  Sgdo  to  minimize  the  geometrical  bias  of  the  discretization.  On  the  basis  of  several  numerical 
simulations  performed  in  this  study,  the  value  dg  =  1.1  seems  to  be  appropriate  to  most  situations.  Finally, 
to  generate  a  statistically  isotropic  random  mesostructure,  the  centers  of  particles  are  placed  throughout 
the  volume  of  the  specimen  one  by  one  (from  the  largest  to  the  smallest)  by  using  a  procedure  introduced 
in  the  concrete  literature  by  Bazant  [5]  and  also  used  by  Cusatis  et  al.  [17].  In  this  procedure,  after 
generating  a  new  particle  position  by  a  random  number  generator,  a  check  is  made  for  possible  overlaps  of 
this  particle  with  previously  placed  particles  and  with  the  surface  nodes.  During  this  phase,  the  surface 
nodes  are  assigned  a  hctitious  diameter  of  Sgdo,  and  a  minimum  distance  of  di/2  +  dj/2  +  (/do  between  the 
centers  of  particles  with  diameters  di  and  dj  is  enforced.  For  (^  =  0  or  very  small  values  of  the  particle 
distribution  tends  not  to  be  statistically  isotropic  and  presents  zones  of  high  particle  density  and  zones 
with  low  particle  density.  On  the  contrary,  for  large  values  of  (,  the  specimen  volume  becomes  saturated 
quickly,  and  particle  placement  cannot  be  completed.  Extensive  numerical  experiments  conducted  by 
the  authors  in  this  research  show  that  (  =  0.2  avoids  volume  saturation  in  most  cases  while  leading  to 
relatively  uniform  particle  distributions.  Figure  Ic  shows  the  particle  system  generated  to  simulate  a 
dog-bone  shaped  specimen. 
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2.3  Geometric  Features  of  the  Interaction  Between  Particles 


The  third  step  of  the  construction  of  concrete  mesostructure  consists  of  dehning  the  topology  of  the 
interaction  among  the  particles.  This  is  obtained  by  Delaunay  tetrahedralization  [21,  2],  which  uses  the 
nodal  coordinates  of  the  particle  centers  as  input  and  gives  a  three-dimensional  mesh  of  tetrahedra  as 
output  that  do  not  overlap,  £11  all  the  volume  of  the  specimen,  and  have  vertices  coinciding  with  the 
given  particle  centers.  In  this  study,  the  Delaunay  tetrahedralization  is  performed  by  using  TetGen  [34]. 
TetGen  implements  a  very  robust  algorithm  for  the  computation  of  conforming  (constrained)  Delaunay 
triangulations  and  allows  modeling  of  non-convex  geometries  such  as  specimens  with  notches  and  cutouts. 
Figure  Id  shows  the  Delaunay  tetrahedralization  of  the  dog-bone  shaped  specimen  whose  particle  system 
is  shown  in  Figure  Ic. 

2.4  Definition  of  Potential  Material  Failure  Locations 

Finally,  the  fourth  step  for  the  characterization  of  concrete  mesostructure  deals  with  the  definition  of 
potential  material  failure  locations  at  the  meso-scale.  As  in  previous  work  by  Gusatis  [17],  it  is  assumed 
here  that  damage/ fracture  initiation  and  evolution  occur  in  the  cement  paste  (or  fine  mortar)  between 
aggregate  pieces,  which  remain  mostly  undamaged  during  the  loading  process. 

In  the  previous  works  by  Gusatis  et  al.  [19,  20],  each  edge  of  the  Delaunay  tetrahedra  was  interpreted 
as  a  connecting  strut  between  two  adjacent  particles,  and  an  effective  area  of  the  strut  was  defined  by 
performing  a  tessellation  of  the  domain  anchored  to  the  Delaunay  tetrahedralization.  This  edge-based 
interaction  is  very  effective  and  allows  a  good  representation  of  concrete  fracturing  behavior  as  well  as 
concrete  failure  under  unconfined  compression.  However,  since  edge-based  interaction  involves  only  the 
kinematics  of  two  particles;  it  does  not  provide  an  accurate  description  of  volumetric  effects  that  are 
important  to  correctly  describe  the  compressive  behavior  of  concrete  under  high  confining  pressure.  As 
will  be  discussed  in  the  next  section,  the  mechanics  of  particle  interaction  in  the  present  study  are  based 
on  an  analysis  of  an  assemblage  of  four  particles  located  at  the  vertices  of  a  tetrahedron.  Gonsequently, 
the  definition  of  the  surfaces  through  which  interaction  forces  are  exchanged  between  particles,  which 
corresponds  to  damage  localization  zones,  is  based  on  the  local  geometry  of  each  tetrahedron  and  the 
corresponding  particles. 


Since  the  tetrahedralization  of  the  particle  centers  is  constructed  through  the  Delaunay  algorithm,  a 
straightforward  tessellation  would  be  the  Voronoi  tessellation.  However,  the  Voronoi  cells  intersect  the 
Delaunay  edges  at  the  mid-points  of  the  edges.  This  feature  is  not  appealing  because,  in  the  general  case 
of  unequal  adjacent  particles,  the  surface  of  the  cells  would  intersect  the  particles  (aggregate  pieces).  More 
in  general,  the  tessellation  of  a  tetrahedron  can  be  obtained  by  a  set  of  triangles  in  which  each  triangle  is 
formed  by  a  point  on  a  tetrahedron  edge  (edge-points,  points  E  in  Figure  le),  a  point  on  a  tetrahedron 
face  (face-points,  points  F  in  Figure  le),  and  a  point  inside  the  tetrahedron  (tet-point,  point  T  in  Figure 
le).  Since  a  tetrahedron  has  six  edges  and  each  edge  is  shared  by  two  faces,  such  a  tessellation  would 
result  in  a  set  of  twelve  triangular  facets. 

The  selection  of  edge-points,  face-points,  and  tet-points  is  somewhat  arbitrary,  but  numerical  experi¬ 
ments  conducted  in  this  study  show  that  the  procedure  outlined  below  tends  to  minimize  the  intersection 
between  the  tessellating  surfaces  and  the  particles.  Such  property  makes  the  representation  of  meso-scale 
crack  paths  consistent  with  the  assumption  that  fracture  initiates  and  propagates  in  the  cement  paste 
and/or  hne  mortar. 

With  reference  to  a  tetrahedron,  (Figure  le)  in  which  vertices  are  labeled  as  1,  2,  3,  and  4,  each  face  is 
labeled  through  the  label  of  the  vertex  opposite  to  that  face,  and  each  edge  is  labeled  through  the  labels 
of  vertices  attached  to  it.  Note  that,  in  Figure  le  as  well  as  in  Figures  If,  Ig,  Ih,  and  2a,  2b,  and  2c,  the 
distance  between  the  particles  was  hctitiously  enlarged  for  illustration  purposes.  Actual  particle  systems, 
such  as  the  one  shown  in  Figure  Ic,  generally  feature  inter-particles  gaps  signihcantly  smaller  than  the 
diameters  of  typical  particles. 

1.  Edge-points  are  dehned  at  midway  of  the  counterpart  of  the  edges  not  belonging  to  the  associated 
particles  (see  point  E12  in  Figure  If  for  the  edge  between  particle  Pi  and  P2). 

2.  Face-points  are  dehned  as  follows.  For  each  face  of  the  tetrahedron,  hrst  three  points  (for  example, 
F41,  F42,  and  F43  for  face  4)  located  on  the  straight  lines  connecting  each  face  vertex  to  to  the 
edge-point  located  on  the  edge  opposite  to  the  particle  under  consideration  are  identihed.  Similarly 
to  the  edge-points,  these  points  are  located  at  midway  of  the  hne  counterpart  not  belonging  to  the 
associated  particles.  In  Figure  Ig,  the  point  F43  associated  with  vertex  P3  and  edge-point  E12  is 
shown.  Then,  the  face-point  is  selected  as  the  centroid  of  these  three  points  (see  points  F  in  Figure 
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le). 


3.  Similarly  to  face  points,  the  tet-point  is  defined  as  the  centroid  T  of  the  points  T^,  T2,  Tg,  and  T4 
identihed  on  the  straight  lines  connecting  each  vertex  of  the  tetrahedron  with  the  face-point  on  the 
face  opposite  to  the  vertex  under  consideration  and  located  at  midway  of  the  line  counterpart  not 
belonging  to  the  associated  particle.  In  Figure  Ih,  the  point  T*4  associated  with  vertex  P4  is  shown. 

By  collecting  all  the  facets  associated  with  one  particle,  one  obtains  a  polyhedral  cell  containing  the 
particle.  Each  cell  has  its  own  random  irregular  shape  (see  Figure  li).  This  property  is  very  important 
to  ensure  a  realistic  representation  of  the  kinematics  of  concrete  mesostructure  and  especially  to  avoid 
excessive  rotations  [36]. 

3  Discrete  Compatibility  and  Equilibrium  Equations 

The  basic  four-particle  tetrahedron,  depicted  in  Figure  Ih,  is  used  here  to  derive  the  governing  equations 
of  the  model.  The  tetrahedron  is  subdivided  into  four  subdomains  Vi  (i  =1,  ...,  4),  each  dehned  by  one 
node  (particle),  a  portion  of  the  three  tetrahedron  edges  attached  to  the  node,  and  the  six  triangular 
tessellation  facets  attached  to  the  three  edges  (see  Figure  2a). 

In  each  subdomain,  the  displacement  held  is  dehned  through  rigid-body  kinematics.  For  x  =  [xi  X2  Xa]"*"  G 
Vi,  one  can  write 


u(x)  =  Ui-F0i  X  (x-Xj)  =  Ai(x)Qi  (4) 

where 


Ai(x) 


1  0  0  0  X3-  X^i  X2i  -  X2 

0  1  0  Xai  -  Xs  0  Xi-  Xu 

0  0  1  X2-X2i  Xu  -  Xi  0 


(5) 


The  vector  Xj  describes  the  position  of  node  i,  and  the  vector  Qj  =  [u]*"  6j]  collects  the  translational, 
uj  =  [uu  U2i  usi],  and  rotational,  6j  =  [9u  92i  6*3i],  degrees  of  freedom  of  node  i  (see  Figure  2a). 

By  using  Equation  4,  it  is  possible  to  dehne  a  displacement  jump  at  the  centroid  C  of  each  facet  in  the 
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tetrahedron 


{uckj  =  Ucj  -  Uc*  (6) 

where  i  and  j  are  the  nodes  adjacent  to  facet  k,  and  ucj  =  and  uci  =  u(x^^)  are  the  displacements 

at  the  facet  centroid  Ck  for  x^^  G  Vj  ,  xj^  G  Vj  (see  Table  1  for  the  permutation  of  indices  i,  j,  and  k). 

By  dividing  the  displacement  jump  by  the  edge  length,  it  is  possible  to  dehne  a  facet  strain  vector 
as  ^|u(7fc],  where  ie  =  ||xj  —  Xj||  =  [(xj  —  Xi)^(xj  —  is  the  length  of  edge  e  (Table  1).  In 

order  to  formulate  a  constitutive  law  featuring  the  classical  tension-compression  asymmetry  of  concrete 
behavior,  the  strain  vector  needs  to  be  decomposed  into  its  normal  and  shear  components.  This 

is  accomplished  with  reference  to  the  projection  of  the  tessellation  facets  into  planes  orthogonal  to  the 
edges  (see  Figure  2b).  The  projected  facets,  as  opposed  as  the  original  ones,  are  used  for  the  dehnition 
of  LDPM  strain  components  to  avoid  non-symmetric  behavior  under  pure  shear.  This  issue  is  clarihed  in 
more  detail  in  Figure  2c  in  which  the  effect  of  a  shear  relative  displacement  between  two  particles,  which 
is  orthogonal  to  the  edge,  is  analyzed.  If  the  facet  is  not  orthogonal  to  the  edge  (see  Figure  2c-left),  the 
shear  relative  displacement  can  cause  either  tension  or  compression  depending  on  its  sign.  Due  to  lack  of 
symmetry  of  concrete  behavior  in  tension  and  compression,  this  leads  to  sign-dependent  meso-scale  shear 
behavior  that,  in  turn,  causes  stress  locking  during  tensile  fracturing  simulations.  Such  a  spurious  effect 
disappears  if  one  formulates  the  mechanical  interaction  between  particles  on  the  projection  of  the  facets 
in  a  plane  orthogonal  to  the  edge  (Figure  2c-right). 

The  strains  components  can  be  then  dehned  as 


^Nk  fj 


(7) 


^Mk  —  f.  — 


(8) 


£Lk  = 


^k  I^Cfcl  _  T:>ikr^ 


(9) 


where  =  (xj  —  Xj)/£e,  and  1^  are  two  mutually  orthogonal  directions  in  the  plane  of  the  projected 
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facets,  =  (l/4)n^Ap(xc-fc),  =  (l/4)m^Ap(xcfc),  and  Bf  =  (1/4)1^ Ap(xcfc),  P  =  i,j-  Equations 

7,  8,  and  9  represent  the  discrete  compatibility  equations  of  the  LDPM  formulation. 

The  meso-scale  constitutive  law,  described  in  the  next  Section,  allows  the  calculation  of  the  normal 
and  shear  stresses  at  each  facet.  Formally,  one  can  write  cr^  =  F{ek,^k)  where  cr^,  and  are  vectors 
collecting  facet  stresses,  strains,  and  internal  variables,  respectively. 

Finally,  the  governing  equations  can  be  completed  by  imposing  the  equilibrium  through  the  Principle 
of  Virtual  Work  (PVW).  The  internal  work  associated  with  a  generic  facet  can  be  expressed  as 

=  ieAkCr'^dSk  =  ^eAkicNkS^Nk  +  C^Mk^^Mk  +  C^Lk^^Lk)  (10) 

where  A^.  is  the  projected  facet  area. 

By  substituting  Equations  7,  8,  9  into  Equation  10,  one  obtains 

5m  =  +  Fj<5Q,  (11) 

F^  =  -^,Ak{aNkB%  +  aMkB%  +  aLuB^t)  (12) 

Fjfc  =  4^fc(o'vfcB^  +  (TMfcB^  +  (TLfcBy^)  (13) 

that  represent  the  nodal  forces  at  node  i  and  j  associated  with  facet  k.  By  summing  up  the  contributions 
of  all  the  facets  and  equating  the  total  internal  work  with  the  total  external  work,  one  can  obtain  the 
discrete  equilibrium  equations  of  the  LDPM  formulation.  It  can  be  shown  that  the  equilibrium  equations 
obtained  through  the  PVW  correspond  exactly  to  the  translational  and  rotational  equilibrium  of  each 
polyhedral  LDPM  cell. 

4  LDPM  Constitutive  Law 

In  this  section,  the  LDPM  constitutive  law,  which  provides  the  relationship  between  the  strain  vector  and 
the  stress  vector  at  the  facet  level,  is  presented.  Hereinafter,  the  subscript  k  indicating  the  generic  facet 
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is  dropped  for  sake  of  simplicity  and  readability  of  the  equations. 


4.1  Elastic  Behavior 

The  elastic  behavior  is  formulated  by  assuming  that  normal  and  shear  stresses  are  proportional  to  the 
corresponding  strains: 


(Ttv  —  aM  —  EtEm'i  o'l  —  E^El  (14) 

where  Epp  =  Eq,  Et  =  aE^,  Eq  =  effective  normal  modulus,  and  a  =  shear- normal  coupling  parameter. 
Eq  and  a  are  assumed  to  be  material  properties  that  can  be  identihed  from  results  of  experimental  tests 
in  the  elastic  regime. 

At  the  macroscopic  level,  concrete  elastic  behavior  is  statistically  homogeneous  and  isotropic.  As  such, 
it  can  be  modeled  effectively  through  the  classical  theory  of  elasticity  characterized  by  Young’s  modulus, 
E,  and  Poisson’s  ratio,  u.  The  relationship  between  the  meso-scale  LDPM  parameters,  a  and  Eo,  and 
the  macroscopic  elastic  parameters,  E  and  u,  can  be  obtained  by  considering  a  limiting  case  in  which  an 
inhnite  number  of  facets  surrounds  one  aggregate  piece.  In  this  case,  the  LDPM  formulation  corresponds 
to  the  kinematically  constrained  formulation  of  the  microplane  model  without  deviatoric/ volumetric  split 
of  the  normal  strain  component  [14].  In  this  case,  one  can  write 


and 


= 


1  -  2z/ 


-E 


2  + 3a 


(15) 


1  -  4z/ 


1  —  q; 


a  = 


V  = 


(16) 


1  +  u  4  +  a 

Figures  3a  and  3b  compare  the  estimates  according  to  Equations  15  and  16  (solid  lines)  to  the  macro¬ 
scopic  average  Poisson’s  ratio  and  Young’s  modulus  obtained  numerically  through  LDPM  simulations 
(circles).  The  LDPM  simulations  were  performed  on  cubic  100-mm-side  specimens  generated  with  the 
following  parameters:  w/c  =  0.5,  a/c  =  6.5,  c  =  300  Kg/m^,  da  =  S  mm,  do  =  4  mm,  and  np  =  0.5.  The 
LDPM  data  points  (circles)  were  obtained  by  setting  Eq  and  a  to  given  values  and  subjecting  the  specimens 
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to  end  displacements  corresponding  to  a  macroscopic  (average)  strain  of  E33  =  2  x  10“^.  The  results  of  the 
simulations  were  post-processed  to  obtain  the  macroscopic  (average)  stress  S33  and  the  macroscopic  (aver¬ 
age)  transverse  strains  En  and  E22  from  which  the  macroscopic  Poisson’s  ratio  and  macroscopic  Young’s 
modulus  were  computed  as  E  =  S33/E33  and  v  =  — 0.5(Eii  -|-  E22)/E33.  The  excellent  agreement  between 
the  numerical  and  the  analytical  results  suggests  that  Equations  15  and  16  can  be  used  conhdently  to 
estimate  the  LDPM  elastic  parameters  from  macroscopic  experimental  measurements  of  Young’s  modulus 
and  Poisson’s  ratio. 

As  it  was  pointed  out  by  Bazant  [6],  Equations  16  do  not  cover  the  entire  range  of  thermodynamically 
acceptable  Poisson’s  ratios  (-1  to  0.5)  since  negative  values  of  a  (and  consequently  negative  values  of  the 
shear  stiffness)  are  associated  with  v  >  0.25.  Although  this  limitation  may  hamper  the  applicability  of 
the  LDPM  framework  to  materials  with  v  >  0.25,  it  does  not  affect  the  modeling  of  concrete  for  which 
Poisson’s  ratio  is  always  about  0.2  or  less.  In  the  microplane  formulation,  the  full  Poisson’s  ratio  range 
can  be  obtained  by  introducing  the  volumetric/deviatoric  decomposition  [14]  of  the  normal  component, 
Sn  =  £v  +  ^Di  and  by  expressing  normal  and  shear  stresses  as  cttv  =  EySv  +  E^iSd,  aM  =  E^iEm, 
ai  =  EjjEl,  where  Ey  =  E / {1  —  2u)  and  Ej:,  =  E/{1  +  u)  are  the  volumetric  and  deviatoric  moduli, 
respectively.  The  same  approach  could  be  used  for  LDPM-type  formulation  [17,  19];  however  in  this  way, 
the  LDPM  capability  of  correctly  simulating  splitting  failure  under  compression  (this  being  one  of  the 
unique  features  of  LDPM)  would  be  completely  lost. 

This  point  can  be  discussed  in  more  detail  considering  Figure  3c,  in  which  a  schematic  representation 
of  the  stress  path  obtained  by  loading  in  compression  a  heterogeneous  material  is  shown.  Heterogeneities 
induce  a  deviation  of  the  compressive  stress  flow  leading  to  the  formation  of  transverse  tensile  stresses  that, 
ultimately  cause  splitting  failure.  In  the  LDPM  formulation  in  which  the  normal  stress  is  proportional  to 
the  normal  strain,  the  transverse  tension  is  automatically  captured.  On  the  contrary,  if  the  deviatoric- 
volumetric  split  is  used,  LDPM-type  formulations  reproduce  the  elastic  solution  of  classical  elasticity  for 
homogeneous  and  isotropic  materials  [30].  In  this  case,  the  stress  compressive  flow  is  undisturbed  from 
the  top  to  the  bottom  of  the  specimen,  and  no  tensile  stresses  can  be  observed.  In  this  situation,  tensile 
failure  is  obviously  not  possible,  and  splitting  cracks  cannot  initiate  and  propagate. 
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4.2  Inelastic  Behavior 

The  LDPM  formulation  of  the  nonlinear  and  inelastic  behavior  aims  at  representing  three  separate  physical 
mechanisms  characterizing  meso-scale  failure  behavior:  1)  fracturing  and  cohesive  behavior  under  tension 
and  tension/shear;  2)  pore  collapse  and  material  compaction  under  high  compressive  stresses;  and  3) 
frictional  behavior. 

4.2.1  Fracturing  Behavior 

Fracturing  behavior  is  characterized  by  tensile  normal  strains  (sat  >  0).  As  previously  developed  by 
Cusatis  [16,  17],  it  is  convenient  to  formulate  fracture  and  damage  evolution  by  a  relationship  between  the 
effective  strain,  e,  and  the  effective  stress,  a,  dehned  as 


By  using  the  effective  strain  and  the  effective  stress,  the  relationship  between  normal  and  shear  stresses 
versus  normal  and  shear  strains  can  be  formulated  through  damage-type  constitutive  equations  (see  deriva¬ 
tion  by  Cusatis  et  ah  [17]): 


S  isf  (XEm 

aN  =  cr — ;  Um  =  ce - ; 

e  e 


ctl  =  cr- 


The  effective  stress  a  is  assumed  to  be  incrementally  elastic,  a  =  EqE,  and  it  is  formulated  to  satisfy 
the  inequality  0  <  a  <  uj),  which  is  enforced  through  a  vertical  (at  constant  strain)  return  algorithm. 
Following  [16,  17],  the  strain-dependent  boundary  abt{£,uj)  can  be  expressed  as 


abt{e,u)  =  ao{u)exp  ^o(^)) 

(Jo(a;) 

in  which  the  brackets  (•)  are  used  in  Macaulay  sense:  (x)  =  maxja;,  0}. 

The  coupling  variable  u  represents  the  degree  of  interaction  between  shear  and  normal  loading  and  is 
dehned  as  [17] 


tana;  = 


En 

y/aST 


(20) 
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where  St  =  total  shear  strain,  and  ar  =  a/ <j\j  +  a\  is  the  total  shear  stress.  The 

boundary  evolves  exponentially  as  a  function  of  the  maximum  effective  strain,  which  is  a  history- 
dependent  variable  dehned  as  Smax  =  \/^Nms,x  +  *^^rmax  where  £v,max(t)  =  max[eAr(r)]  and  £r,max(t)  = 
max[£7’(r)]  are  the  maximum  normal  and  total  shear  strains,  respectively,  attained  during  the  loading 

T<t 

history.  It  is  worth  noting  that  in  absence  of  unloading,  e^ax  =  £• 

The  function  ao{u)  is  the  strength  limit  for  the  effective  stress  and  is  dehned  as 

—  sin(a;)  -|-  A/sin^(a;)  +  Aacos^{u) 

^o(^)  ^  2  {  \  I  2 

2acos^[ijj)/rg^ 

in  which  =  ag/ut  is  the  ratio  between  the  shear  strength  (cohesion),  Us,  and  the  tensile  strength,  at-  In 
the  stress  space  a^  —  aT,  Equation  21  is  a  parabola  with  its  axis  of  symmetry  coinciding  with  the  ciAr-axis 
(solid  curve  in  Figure  4a). 

The  exponential  decay  of  the  boundary  au  starts  when  the  maximum  effective  strain  reaches  its  elastic 
limit  £o(^)  =  <^o(i^)/-Eo,  and  the  decay  rate  is  governed  by  the  post-peak  slope  (softening  modulus),  which 
is  assumed  to  be  a  power  function  of  the  internal  variable  cu: 


Ho{u)  =  Ht  ‘  (22) 

Equation  22  provides  a  smooth  transition  from  softening  behavior  under  pure  tensile  stress,  Hq{'k/2)  = 
Ht,  to  perfectly  plastic  behavior  under  pure  shear,  Hq{0)  =  0.  In  order  to  preserve  the  correct  energy 
dissipation  during  meso-scale  damage  localization  [4],  the  softening  modulus  in  pure  tension  is  expressed 
as  Ht  =  2Eq/  [it!^  —  1),  where  Gt  is  the  meso-scale  fracture  energy,  it  =  2EttGtla\,  and  I  is  the  length  of 
the  tetrahedron  edge  associated  with  the  current  facet. 

The  dashed  curve  in  Figure  4a  represents  the  strength  domain  (computed  from  the  boundary  a^) 
corresponding  to  £max  =  4(Ti/Eo.  As  one  can  see,  as  damage  progresses,  the  strength  domain  shrinks 
and  tends  to  become  concave.  This  characteristic  is  important  in  order  to  correctly  simulate  the  different 
ductilities  in  tension  and  compression  at  the  macro-scale.  It  is  also  worth  noting  that,  in  damage  mechanics 
(contrary  to  plasticity),  the  lack  of  convexity  of  the  strength  domain  does  not  violate  thermodynamics 
restrictions  on  the  energy  dissipation,  and  as  such,  it  is  a  completely  acceptable  feature  of  the  formulation. 

Figure  4b  shows  the  normal  and  shear  stress  versus  strain  relationships  for  a;  =  0  (pure  shear),  a;  =  7r/2 
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(pure  tension),  and  u  =  tt/S.  The  plots  were  generated  by  setting,  for  illustration  purposes,  £=  10  mm, 
a  =  0.25,  Ef)  =  30000  MPa  ,  at  =  3  MPa,  ag  =  4.5  MPa,  it  =  100  mm,  and  Ut  =  0.2. 

Finally,  the  fracturing  formulation  needs  to  be  completed  by  unloading-reloading  rules  to  simulate 
cycling  loadings.  Figure  4c  illustrates  the  unloading- reloading  rule  adopted  in  this  work  in  terms  of 
effective  stress  versus  effective  strain  relationship.  Let’s  assume  that  unloading  occurs  after  the  effective 
strain  increased  continuously  from  zero  to  a  certain  value  ^max-  The  effective  stress  decreases  elastically 
until  it  reaches  a  zero  value  and  remains  constant  at  zero  for  further  decreases  of  the  effective  strain. 
During  reloading,  the  effective  stress  remains  zero  until  the  effective  strain  reaches  the  reloading  strain 
limit,  Etr,  and,  beyond  this  point,  the  behavior  is  incrementally  elastic.  The  reloading  strain  limit  is  dehned 
as  Etr  =  kt{Ema.x  —  o'bt/Eo),  where  kt  is  assumed  to  be  a  material  parameter.  The  parameter  kt  governs  the 
size  of  hysteresis  cycles  and,  consequently,  the  amonnt  of  energy  that  the  material  can  dissipate  during 
cycling  loading.  For  kt  =  1,  the  dissipated  energy  is  zero,  whereas  for  kt  =  0,  it  is  a  maximnm  and  eqnal 
to  AiahtEm^x  for  a  given  E^ax  on  a  facet  of  projected  area  A  and  associated  with  a  lattice  edge  length  i. 

4.2.2  Pore  Collapse  and  Material  Compaction 

Under  high  compressive  hydrostatic  deformations,  concrete  behavior  exhibits  strain-hardening  plasticity. 
The  plastic  behavior  is  characterized  by  an  initial  phase  in  which  micro-scale  and  meso-scale  pores  collapse 
nnder  load,  and  a  later  phase,  in  which  the  walls  of  completely  collapsed  pores  are  in  contact  leading  to  a 
signihcant  densihcation  of  the  material  [25].  In  terms  of  stress-strain  response,  the  hrst  phase  is  associated 
with  a  sndden  decrease  of  the  stiffness  that  is  later  regained  in  the  second  phase  (rehardening).  Experiments 
show  [3]  that  after  the  densihcation  phase,  both  the  tangent  plastic  stiffness  and  the  nnloading  elastic 
stiffness  are,  in  some  cases,  even  larger  than  the  initial  elastic  stiffness.  Finally,  rehardening  is  limited  (or 
even  negligible)  in  the  presence  of  signihcant  deviatoric  deformations,  which  typically  leads  to  a  horizontal 
plateau  in  the  measnred  macroscopic  stress  versus  strain  cnrves. 

LDPM  constitutive  law  simulates  the  phenomena  discussed  above  through  a  strain-dependent  normal 
bonndary  limiting  the  compressive  normal  stress  component  at  the  facet  level.  This  compressive  boundary, 
abc{£Di£v)i  is  assnmed  to  be  a  function  of  the  volumetric  strain  Ey  and  the  deviatoric  strain  ed-  The 
volnmetric  strain  is  computed  at  the  tetrahedron  level  as  Ey  =  iV  —  Vq) /Vq  where  V  and  Vq  are  the  cnrrent 
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and  initial  volume  of  the  tetrahedron,  respectively.  The  current  volume  is  computed  by  neglecting  the  effect 
of  nodal  rotations.  A  more  rigorous  definition  of  volumetric  strain  in  the  framework  of  discrete  models  and 
accounting  particle  rotations  is  provided  by  Cusatis  and  Schauffert  [?].  In  each  LDPM  tetrahedron,  all 
twelve  facets  are  assumed  to  be  subject  to  the  same  volumetric  strain  whereas  each  facet  is  characterized 
by  a  different  value  of  the  deviatoric  strain  calculated  by  subtracting  the  volumetric  strain  from  the  normal 
strain:  ed  =  £n  —  £v-  The  definition  of  the  volumetric  and  deviatoric  strains  are  completely  equivalent  to 
the  same  quantities  defined  at  each  microplane  in  the  microplane  model  formulation  [6]. 

For  a  constant  deviatoric-to- volumetric  strain  ratio,  r^v  =  Ed/sv,  the  compressive  boundary  is  as¬ 
sumed  to  have  an  initial  linear  evolution  (modeling  pore  collapse  and  yielding)  followed  by  an  exponential 
evolution  (modeling  compaction  and  rehardening).  One  can  write 


[  CTco  for  -Ev  <0 


O'bci^D,  £v)  —  S 


o'co  +  {—£v  —  £c'd)Hc{r-ov) 


for  0  <  —Ev  <  £ci 


o-ciirDv)  exp  [{-Ev  -  £ci)Hc{rDv) / c^dirDv)]  otherwise 


(23) 


where  aco  is  the  meso-scale  yielding  compressive  stress,  Eco  =  ctcoZ-Eq  is  the  compaction  strain  at  the  onset 
of  pore  collapse,  Hc{rjyv)  is  the  initial  hardening  modulus,  Ed  =  ^co^co  is  the  compaction  strain  at  which 
rehardening  begins,  Kco  is  the  material  parameter  governing  the  onset  of  rehardening,  and  Udi'^Dv)  = 
o'co  +  (^ci  —  £co)Hc{rDv)- 

For  increasing  rjjv,  the  slope  of  the  initial  hardening  modulus  needs  to  tend  to  zero  in  order  to  simulate 
the  observed  horizontal  plateau  featured  by  typical  experimental  data.  This  can  be  achieved  by  setting 


Hc{rDv) 


1  +  «c2  {tdV  -  «cl) 


(24) 


where  FfcO;  ^ci,  and  Kc2  are  assumed  to  be  material  parameters.  For  compressive  loading  {ejv  <  0),  the 
normal  stress  is  computed  by  imposing  the  inequality  —(Jbc{£D,£v)  <  o'v  <  0.  Within  the  boundaries  of 
this  inequality,  the  behavior  is  assumed  to  be  incrementally  elastic:  &n-  =  E]Vc£n-  In  order  to  model  the 
increased  stiffness  during  unloading,  the  loading-unloading  stiffness  E^c  is  defined  as 
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Eq  for  -aN  <  CTco 


(25) 


Enc  =  \ 

Ed  otherwise 

where  Ed  is  the  densihed  normal  modulus.  For  loading  processes  at  constant  rov,  we  have  ey  =  + 

r£)y).  In  this  case,  the  boundary  in  Equation  23  can  be  expressed  as  a  function  of  the  normal  strain  ejv, 
shown  in  Figure  4d,  where  the  solid  and  dashed  curves  represent,  respectively,  the  compressive  normal 
stress  versus  strain  relationship  for  r^y  =  0  and  r^y  =  1.1.  For  the  case  of  roy  =  0,  the  unloading¬ 
reloading  rule  is  also  shown.  For  illustration  purposes,  we  assumed  parameter  values  of  Eq  =  60000  MPa, 
(Tco  =  100  MPa,  Hco/Eq  =  0.6,  Kco  =  ScilSco  =  4,  Ed/E^  =  2,  Kci  =  1,  and  k,c2  =  5. 

4.2.3  Frictional  Behavior 

Finally,  in  the  presence  of  compressive  stresses,  the  shear  strength  increases  due  to  frictional  effects.  As 
often  done  in  the  literature,  frictional  phenomena  can  be  simulated  effectively  through  classical  incremental 
plasticity  [22].  Incremental  shear  stresses  can  be  calculated  as  au  =  Et^Em 

where  the  plastic  strain  increments  are  assumed  to  obey  the  normality  rule  =  \dip/daM  and  = 
Xdip/daL- 

The  plastic  potential  can  be  expressed  as  99  =  ~  (^bs{(EN)  in  which  the  shear  strength  ahs  is 

formulated  with  a  nonlinear  frictional  law: 

o'bsio'N)  =  0's  +  (/^o  —  f^oo)o'm  —  ^^ooO^N  —  (/xq  ”  f^oo)o'm  exp  {on/ono)  (26) 

where  (jg  is  the  cohesion  (previously  introduced  in  this  paper),  /tq  and  /ioo  are  the  initial  and  hnal  internal 
friction  coefficients,  respectively,  and  a  no  =  is  the  normal  stress  at  which  the  internal  friction  coefficient 
transitions  from  /tq  to  /ioo-  As  shown  in  Figure  4e,  the  frictional  law  in  Equation  26  tends  to  the  asymptote 
o^s  +  (ho  ~  1^oo)ono  ~  I^ooOn  characterized  by  a  slope  equal  to  /ioo  and  whose  intercept  with  the  a-r-axis 
increases  with  increasing  ono-  The  parameter  aNo  basically  governs  the  extent  of  the  nonlinearity  of  the 
shear  boundary.  The  classical  linear  (Coulomb-type)  frictional  law  with  slope  /io  or  /ioo  is  obtained  by 
setting  a  NO  =  00  or  (Tato  =  0,  respectively.  Finally,  equations  governing  the  shear  stress  evolution  must 
be  completed  by  the  loading-unloading  conditions,  which  can  be  expressed  through  the  classical  Kuhn- 
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Tucker’s  conditions:  (p\  <  0  and  A  >  0.  Typical  shear  versus  strain  relationship  is  shown  in  Figure 
4f. 

5  Numerical  Implementation  and  Stability  Analysis 

The  LDPM  computational  framework  formulated  in  this  paper  was  implemented  into  the  MARS  code 
[27],  which  is  a  multi-purpose  structural  analysis  program  based  on  a  modern  object-oriented  architecture. 
MARS  performs  structural  analysis  by  an  explicit  dynamic  algorithm  (based  on  a  central  difference  scheme) 
and  is  very  effective  in  the  management  of  the  various  computational  entities  (nodes,  hnite  elements, 
loads,  etc.)  making  possible  the  numerical  simulation  of  very  large  systems  even  on  regular  desktop 
computers.  These  features  are  particularly  attractive  for  the  research  presented  in  this  paper  since  LDPM 
calculations  often  involve  several  thousands  of  degrees  of  freedom.  In  addition,  the  explicit  character  of 
the  computational  scheme  implemented  in  MARS  makes  it  very  advantageous  because  is  not  affected  by 
the  convergence  problems  that  implicit  schemes  often  have  in  handling  softening  behavior. 

Explicit  algorithms,  however,  are  not  unconditionally  stable  and  require  an  accurate  evaluation  of 
the  numerical  stability  of  the  numerical  simulations.  In  the  elastic  regime,  the  stability  condition  can 
be  expressed  as  At  <  2/a;max  where  cUmax  represents  the  highest  natural  frequency  of  the  computational 
system.  It  can  be  shown  [8]  that  cUmax  <  max(a;7)  where  uj  are  the  natural  frequencies  of  the  individual 
unrestrained  elements  composing  the  mesh  used  in  the  simulation.  Based  on  this  observation,  the  stable 
time  step  for  the  LDPM  can  be  estimated  by  computing  the  natural  frequencies  of  each  LDPM  tetrahedra. 
This  requires  solving  the  eigenvalue  problem  det(K  —  ca^M)  =  0  where  K  is  the  stiffness  matrix  and  M  is 
the  mass-matrix. 

The  elastic  energy  associated  with  a  generic  facet  k  is 

Ut  =  +  Ereh  +  Erel,,)  =  (K‘  +  K‘.)Qi  +  (Kf.  +  KpQ.  (27) 

where  +  Et(b£)TB1?  +  )TBf  (p,  q  =  *, j). 

The  overall  tetrahedron  stiffness  matrix  is  obtained  by  assembling  strain  energy  contributions  of  all 
twelve  facets: 
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(28) 
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12 


k=l 


k=l 


K?,- 

Kf,- 

K^. 

where  the  symbol  'Yh  is  used  to  identify  the  assemblage  operation.  The  kinetic  energy  associated  with  a 
generic  facet  k  can  be  subdivided  into  two  terms  associated  with  the  two  nodes  {i  and  j,  see  Table  1) 
adjacent  to  the  facet:  Tk  =  Tki  +  %j.  Each  individual  term  reads 


7kp 


'Vi 


kp 
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(30) 


and  p  =  c{w /c  +  a/c)  +  PairVair  is  material  density,  pair  is  1.2  Kg/m^  (at  sea  level  and  at  15°C),  and  Vkp 
is  the  volume  identihed  by  the  facet  k  and  the  node  p.  The  various  terms  appearing  in  the  matrix 
are  hrst  order  (symbols  S)  and  second  order  (symbols  X)  moments  of  the  volume  Vkp  about  the  axes  of  a 
cartesian  system  of  reference  with  origin  at  node  p. 

Similar  to  the  stiffness  matrix,  the  overall  tetrahedron  mass  matrix  is  obtained  by  assembling  the 
contributions  of  the  twelve  facets 
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12 


M  =  ^  ^ 


fc=i 


k=l 


Mf  0 
0 


(31) 


The  mass  matrix  in  Equation  31  is,  in  general,  not  diagonal.  In  order  to  take  full  advantage  of  the 


explicit  algorithm  used  in  this  study,  a  diagonalized  version  of  M  is  obtained  by  simply  discarding  the 
non-diagonal  terms.  It  is  worth  noting  that  by  using  the  diagonalized  matrix,  the  mass  term  of  a  certain 
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particle  node  is  exact  if  the  particle  node  is  the  mass  centroid  of  the  LDPM  cell.  This  is  approximately 
the  case  for  many  particle  nodes  in  the  interior  of  typical  LDPM  systems. 

6  Concluding  Remark 

The  Lattice  Discrete  Particle  Model  (LDPM)  formulated  in  this  paper  has  several  unique  features  and 
potential  capabilities,  which,  however,  need  to  be  verified  by  performing  numerical  simulations  and  com¬ 
paring  numerical  results  with  experimental  data.  Part  II  of  this  study,  which  follows,  presents  an  extensive 
calibration  and  validation  of  LDPM. 
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facet 
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e 

node 

i 

node 
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1,2 
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2 

3,4 

2 

1 

3 

5,  6 

3 

2 

3 

oo 

4 

2 

4 

9,  10 

5 

3 

4 

11,  12 

6 

1 

4 

Table  1:  Facet,  edge,  and  node  indices. 
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Figure  1:  a)  Probability  distribution  function  for  particle  size  generation;  b)  Theoretical  (solid  curve)  and 
numerical  (circles)  sieve  curve;  c)  Particle  system  for  a  typical  dogbone  specimen;  d)  Tetrahedralization 
for  a  typical  dogbone  specimen;  e)  Tessellation  of  a  typical  LDPM  tetrahedron  connecting  four  adjacent 
particles;  f)  Edge-point  dehnition;  g)  Face-point  dehnition;  h)  Tet-point  dehnition;  i)  LDPM  cells  for  two 
adjacent  aggregate  particle. 
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Figure  2: 
facets;  c) 
projected 


a)  Counterpart  of  a  tetrahedron  associated  to  one  particle;  b)  Original  and  projected 
Effect  of  meso-scale  pure  shear  loading  on  one  LDPM  original  facet  (left)  and  one 
facet  (right). 


LDPM 

LDPM 


Figure  3:  a)  Relationship  between  macroscopic  Young’s  modulus  and  LDPM  elastic  parameters  ([-]  means 
that  the  correspondent  quantity  is  adimensional) ;  b)  relationship  between  macroscopic  Poisson’s  ratio  and 
LDPM  coupling  parameter;  c)  Idealization  of  the  effect  of  compression  on  heterogeneous  materials. 
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Figure  4:  LDPM  constitutive  laws,  a)  Virgin  and  damaged  shear  strength  as  function  of  positive  normal 
stresses;  b)  Typical  stress  versus  strain  curves  at  the  LDPM  facet  level;  c)  Unloading-reloading  model; 
d)  Typical  normal  stress  versus  normal  strain  curves  in  compression;  e)  Shear  strength  as  a  function  of 
negative  normal  stresses  (frictional  behavior);  f)  Typical  shear  stress  versus  shear  strain  curve. 
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